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1 Introduction 

The study of self-gravitating stellar systems has provided important hints to 
develop tools of analytical mechanics. We may cite the ideas of Jeans pQ about 
the relevance of conserved quantities in describing the phase-space structure of 
large N-body systems and his introduction of the concept of isolating integral. 
Later important contributions are those of Contopoulos [3J, who applied a 
direct approach to compute approximate forms of the isolating integrals of 
motion, of Henon & Heiles [8] with a paradigmatic example of non-integrable 
system derived from a simple galactic model and of Hori [7] , who introduced 
the theory of Lie transforms in the field of canonical perturbation theory. 
These and other cues contributed to the body of methods and techniques 
that we use today to study regular and chaotic rdynamics of non-integrable 
systems. 

The direct approach applied by Contopoulos aims at solving the equation 
for the conserved quantity along the classical procedure developed early in the 
last century [3] . The method of the Lie transform 7] , subsequently improved 
by several authors [HI E3 El EH EU > has some technical advantages and has 
gradually become a standard method in the perturbation theory of Hamil- 
tonian dynamical systems [4]. However, its application in galactic dynamics, 
with a few remarkable exceptions Ell El] , has not been as systematic and 
productive as it could be. 

Hamiltonian normal forms constructed in this way [161 1171 118j are a pow- 
erful tool to investigate the orbit structure of galactic potentials and to gather 
several qualitative informations concerning the near integrable dynamics be- 
low the stochasticity threshold (if any) of the system. Results obtained in the 
same class of systems by the averaging method [19] are easily overtaken. As a 
matter of fact, with a normal form truncated to an order sufficient to incorpo- 
rate the main resonance, one can also make reliable quantitative predictions. 
In the present contribution we review how to exploit detuned resonant normal 
forms to extract information on several aspects of the dynamics in systems 
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with self-similar elliptical equipotentials. In particular, using energy and el- 
lipticity as parameters, we compute the instability thresholds of axial orbits, 
bifurcation values of low-order boxlets and phase-space fractions pertaining 
to the families around them. We also show how to infer something about the 
singular limit of the potential. 

A remarkable side-effect of expressing the stability-instability threshold 
as a series expansion, is that its predictive ability goes well beyond the radius 
of convergence of the perturbing expansion. Exploiting asymptotic properties 
of the series constructed via the normal form [201 HI] > we mav try to estimate 
an optimal truncation order. 



2 The Hamiltonian normal form 

The subject of our investigation is the class of 2-dof natural systems 

H (P, r) = \<J>1+ Pl/q) + V(s(x, y)). (1) 
V is a uniformly increasing function of the variable 

s = x 2 +y 2 /q, (2) 

with an absolute regular minimum (1^(0, 0) = ^'(0, 0) =0), so that the energy 
E may take any non-negative value. Two simple examples are 

V L = Uog{l + s ), (3) 
V c = VTT~s-l. (4) 

The parameter q gives the "ellipticity" of the equipotentials and ranges in the 
interval 

0.6<g<l. (5) 

Lower values of q can in principle be considered but correspond to an unphys- 
ical density distribution if V is a gravitational potential. Values greater than 
unity are included in the treatment by reversing the role of the coordinate 
axes. With respect to the standard 'physical' notation, the scaling transfor- 
mation 

Py — > Vqpy, v — > y/V? ( 6 ) 

is implicit in the Hamiltonian written in the form 
2.1 Series expansions 

To investigate the dynamics of system ((TJ), we look for a new Hamiltonian 
given by the series expansion in the new canonical variables P, R, 
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K(P,R) = J2K n (P,R), (7) 

n=0 

with the prescription that 

{H o ,K} = 0. (8) 

In these and subsequent formulas we adopt the convention of labeling the first 
term in the expansion with the index zero: in general, the 'zero order' terms 
are quadratic homogeneous polynomials and terms of order n are polynomials 
of degree n + 2. The zero order (unperturbed) Hamiltonian, 

H (P, R)^K = 1{P% + X 2 ) + j-(P 2 + F 2 ), (9) 

with unperturbed frequencies u)\ = 1 and 0J2 = l/q, is expressed in terms of 
the new variables found at each step of the normalizing transformation. It is 
customary to refer to the scries constructed in this way as a "Birkhoff" normal 
form [5]. The presence of terms with small denominators in the expansion, 
forbids in general its convergence. It is therefore more effective to work since 
the start with a resonant normal form [6], which is still non-convergent, but 
has the advantage of avoiding the small divisors associated to a particular 
resonance. To catch the main features of the orbital structure, we therefore 
approximate the frequencies with a rational number plus a small "detuning" 

Ul rri\ 

— = q = h 0. (10) 

We speak of a detuned (1711/1712) resonance, with mi + m-i the order of the 
resonance. In order to implement the normalization algorithm, also the orig- 
inal Hamiltonian JTJ) has to be expressed as a series expansion around the 
equilibrium: performing the rescaling 

n . = rn2H = 

we redefine the Hamiltonian as the series 

n = Y, Hk =2 [mi ^ + x2 ) + m 2 (Pi + V 2 )] + irn 2 S(p 2 x + x 2 ) + h(q)s k+ \ 

k=0 k=l 

(12) 

with expansion coefficients bk depending only on the ellipticity in view of the 
restriction imposed by the choice of the potentials. The procedure is now that 
of an ordinary resonant "Birkhoff-Gustavson" normalization [22J, [23] with 
two variants: the coordinate transformations are performed through the Lie 
transform and the detuning quadratic term is treated as a term of higher order 
and put into the perturbation. This is analogous to the strategy of the 'nearly 
resonant construction' of Contopoulos and Moutsoulas [24] in the context of 
the direct approach and is implemented in the program by Giorgilli [25] , 
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2.2 Lie transform normalization 

Considering a generating function g, the new coordinates P, R result from 
the canonical transformation 

(P,R) = M g (p,r). (13) 
The Lie transform operator M g is defined by [3] 

oo 

M g = Y,M k (14) 

fc=0 

where 

fc . 

M = 1, M k =J2 J.I ■, U/ ;• (15) 

3=1 

The functions gj are the terms in the expansion of the generating function 
(go = 1) and the linear differential operator L g is defined through the Poisson 
bracket, L g (-) = {g, ■}. 

The terms in the hew Hamiltonian are determined through the recursive 
set of linear partial differential equations 0] 

n-1 

K n =H n + Y J M n-j'H.j, 7i=l,2,... (16) 

3=0 

'Solving' the equation at the n-th step consists of a twofold task: to find 
K n and g n . We observe that, in view of the reflection symmetries of the 
Hamiltonian (|TJ) , the chain l|16p is composed only of members with even index 
and so the normal form itself is composed of even-index terms only. The 
unperturbed part of the Hamiltonian, Tlo, determines the specific form of the 
transformation. In fact, the new Hamiltonian K is said to be in normal form 
if, analogously to (|8|) , 

{Ho,K} = 0, (17) 

is satisfied. The function 

T = K - Ho (18) 

can be used as a second integral of motion. For practical applications (for 
example to compare results with numerical computations) it is useful to ex- 
press approximating functions in the original physical coordinates. Inverting 
the coordinate transformation, the new integral of motion can be expressed 
in terms of the original variables. Denoting it as the power series 

oo 
n=0 

its terms can be recovered by means of the set of equations 
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^ M n - 3 [Hj - I,] -K n = 0, n > 1, (20) 

that is obtained from l[T6|) and (|18jl by exploiting the nice properties of the 
Lie transform with respect to inversions [I]. 

2.3 Effective order of the detuning 

We have to discuss how to treat the detuning term: it is considered as a higher- 
order term and the most natural choice is to put it into H.2- However, there is 
no strict rule for this and one may ask which is the most 'useful' choice, always 
considering that applications are based on series expansions with coefficients 
depending on q. We remark that, different choices of the effective order, say 
d, of the detuning, lead to different terms of higher order in the normal form. 
We also observe that, whatever the choice made, the algorithm devised to 
treat, step by step, the system ([To]) must be suitably adapted to manage with 
polynomials of several different orders. In practice, since at each step the 
actual order of terms associated to detuning is lower than the corresponding- 
effective order, the algorithm is adapted by incorporating routines already 
used at previous steps. 

In practice, at step say j, we have an equation of the form 

Kj = Hj + Aj + 8Bj_ d + 8 2 B^ 2d + ... + L 9] (H ), (21) 

where A4 , Bi are homogeneous polynomials of degree i + 2 coming from pre- 
vious steps. As usual, the algorithm is designed to identify in all terms with 
the exclusion of 

L gj (n )=-L no ( gj ), (22) 

monomials in the kernel of the linear operator L-j-i . These monomials are used 
to construct Kj : the remaining terms are used to find gj in the standard way. 
It is clear that both the normal form and the generating function are affected 
by the effective order of the detuning term. 

In both cases (|3l4p investigated |2S] , with the detuning treated as a term 
of order 2, the next appearance of a related term is in K e . Rather, if it is 
treated as a term of order 4, the next appearance of a related term is in Kg. 
Truncating at order 6 (polynomials of degree 8) is therefore sufficient to make 
a comparison with other predictions not sensitive to the detuning. 

2.4 Structure of the normal form 

In principle, the recursion process to solve the system (|16l) can be carried out 
to arbitrary order. In practice we have to truncate it at some finite order TV. 
The ideal choice would be an optimal truncation order N opt , to evaluate which 
one can work with the formal integral (|I9[) : minimizing its failure to commute 
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with the Hamiltonian, one truncates the series at the order giving the best 
conservation [J_21 IHH 121] • However, in general this is a costly procedure. 

On the other hand, a very conservative strategy can be that of truncating 
at the lowest order adequate to convey some non-trivial information on the 
system. In the resonant case, it can be shown [27] that the lowest order to be 
included in the normal form in order to capture the main effects of the m\jm 2 
resonance with double reflection symmetries is iV m i n — 2 x (mi + 11J2 — I). 
Truncating at this level is enough to study the resonance and the main periodic 
orbits associated to it [TBI 02]. Using ' action-angle '-like variables J, 6 defined 
through the transformation 

X = V^jicosfli, P x = v/2Ji~sm<9i, (23) 
Y = y/2J^cos6 2 , P y = V2J2sin6> 2 , (24) 

the typical structure of the doubly-symmetric resonant normal form truncated 
at iV min is [6] 

mi+m 2 

K = m 1 J 1 +m 2 J 2 + V {k) {Jx, J 2 )+aJr 2 J 2 mi cos[2(™ 2 6' 1 -m 1 6l 2 )], (25) 

k=2 

where V {k) are homogeneous polynomials of degree k whose coefficients may 
depend on S and a is a constant. In these variables the second integral is 

£ = mi Ji + m 2 J 2 (26) 

and the angles appear only through the resonant combination 

ij) = m 2 6 1 - mi 6 2 . (27) 

Introducing its conjugate variable 

1Z = m 2 Ji — m± J2, (28) 

the new Hamiltonian can be expressed in the reduced form K(1Z, tp; £, S) that 
is a family of 1-dof systems parametrized by £ and S. 

3 Applications 

We now analyze a sample of problems that can be addressed and solved with 
the tools developed so far. The main periodic orbits and the families of quasi- 
periodic orbits parented by them give the structure of the phase space, at least 
in the regular regime. The study of existence and stability of normal modes 
and periodic orbits in general position admitted by ([25]) or its higher-order 
generalization is therefore of outmost importance. 
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3.1 Stability of normal modes 

In systems of the form ([I} the orbits along the symmetry axes are simple 
periodic orbits. It can be readily verified that these orbits correspond to the 
two solutions in which either J\ or Ji vanish. If the axial orbit is stable it 
parents a family of 'box' orbits. A case that is both representative of the state 
of affairs and useful in galactic applications is that of the stability of the x- 
axis periodic orbit (the 'major-axis orbit', if q is in the range l|5|)). Among 
possible bifurcations from it, the most prominent is usually that due to the 
1:2 resonance between the frequency of oscillation along the orbit and that 
of a normal perturbation, producing the 'banana' and 'anti-banana' orbits 
[28] . The inclusion of detuning allows one to catch the passage of the system 
through the resonance due to the nonlinear coupling between the two degrees 
of freedom: the strength of the coupling depends on energy and we expect 
that the onset of the resonance is described by one (or more) curves on the 
(S, i?)-plane. To investigate this problem in the potentials (|3I4|) , we construct 
the normal form with mi = 1,TO2 = 2 and study the nature of the critical 
points of the function K^> = K + fJ-Ti-a, where \i has to be considered as a 
Lagrange multiplier to take into account that there is the constraint 7io = £ ■ 
The condition for a change in the nature of the critical point corresponding 
to the normal mode is given by the solutions of the algebraic equation 

dct [d 2 !^ (£,<$)] |. /2=0 = (29) 

of degree N in £ : each transition of the kind extremum — ► saddle is equivalent 
to the onset of an instability and to the bifurcation of the banana (or of the 
anti-banana) . 

However, in order to get a form usable in comparison with other results (for 
example coming from a numerical treatment) it is necessary to use a 'physical' 
energy variable rather than the parameter £. The conversion is possible if the 
physical energy E appears explicitly 17J. According to the rescaling (fTTj) , we 
assume that ra-iqE is the constant 'energy' value assumed by the truncated 
Hamiltonian K. In the present instance ni2 = 2 so that, on the x-axis orbit, 
the new Hamiltonian is a series of the form 

K = 2q£ + cq£ 2 + ... = 2qE. (30) 

This series can be inverted to give 

£ = E-^E 2 + ... (31) 

and this can be used in the treatment of stability to replace £ with E. Recalling 
that, in this case, (fit))) gives q = 1/2 + S, every solution can be expanded as 



N/2 

E alt (6) = J2c k 6 k 
fc=i 



(32) 
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Table 1. Coefficients in the expansion (|32|l with N = 14 for the logarithmic poten- 
tial (banana, 2nd column and anti-banana, 3rd column) and the conical potential 
(banana, 4th column and anti-banana, 5th column). 





Potential Vl 


Potential Vc 


k 


Banana 


Anti — banana 


Banana 


Anti — banana 


1 


8 


8 


16 


16 


2 


20 
3 


28 
3 


248 
3 


536 
3 


3 


268 
9 


460 
9 


3608 
9 


18584 
9 


4 


1724 
27 


3928 
27 


43328 
27 


657848 
27 


5 


79184 
405 


267404 
405 


525704 
81 


23668304 
81 


6 


567178 
1215 


510200857 
405 


28118794 
1215 


4304374384 
1215 


7 


30991946 
25515 


615376795556 
8505 


309430864 
3645 


31575390356 
729 



and in this form they can be used for quantitative predictions. 

In Table [TJ we list the coefficients of the series ([32j) giving these bifurca- 
tions for the logarithmic potential ([3|) and the 'conical' potential (j4|). They 
have been obtained [26) with a normal form truncated at order N — 14 and 
with the detuning treated as a term of order 2. There is a complete agree- 
ment with the analytical approach based on the Poincare-Lindstedt method 
|29j and, as discussed below, there is a striking agreement with the numeri- 
cal approach based on the Floquet method. The agreement of all fractional 
coefficients is complete up to N/2 — 7. On the other hand, if the detuning 
is treated as a term of order d = 4 or greater, we get a disagreement in the 
coefficients starting from C3 . This result confirms the analysis made above on 
the 'propagation' of the detuned terms in the normal form and show that the 
choice d = 2 is the optimal one. 

What is remarkable in the quality of the prediction with regard to 'exper- 
imental' numerical data is that numerical computations are performed with 
the exact logarithmic (or conical) potentials (|3I4[) , whereas the analytical pre- 
dictions are, in any case, based on the series expansions of these potentials 
that appear in (| 12(1 with limited convergence radii. The reliability of these 
predictions in a range wider than foreseen can be explained if we interpret the 
series of the form (|32p as asymptotic series and evaluate their truncations by 
computing the successive partial sums 

n 

E n (q)=Y. Ck5k ' n = l,...,N/2. (33) 

k=l 

Minimizing the difference between the 'exact' value and its approximations 
provides an estimate of the optimal truncation. As an example, in Table O 
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Table 2. Subsequent truncations of expansion (|32|1 with N = 14 for the logarithmic 
potential (banana). Eb is the value obtained by means of the Floquet method. 





8 


n 


0.1 


0.2 


0.3 


0.4 


1 


n snnnnn 


i rsnnnn 




O.iUUUU 


2 


0.733333 


1.33333 


1.80000 


2.13333 


3 


0.763111 


1.57156 


2.60400 


4.03911 


4 


0.756726 


1.46939 


2.08680 




5 


0.758681 


1.53196 


2.56190 




6 


0.758214 


1.50208 


2.22160 




7 


0.758336 


1.51763 


2.48724 




E B 


0.758 


1.513 


2.401 


3.646 



we report these partial sums for the banana bifurcation in the logarithmic 
potential ([3]), with 0.1 < S < 0.4 (0.6 < q < 0.9) and compare them with 
the values obtained by means of the Floquet method [28l [17] given in the last 
row. The numerical values of the partial sums are given with 6 digits just to 
show more clearly the asymptotic behaviour: we can see that, up to S = 0.3, 
the predictions are apparently still (slowly) converging at n = 7. Only at the 
rather extreme value 5 = 0.4 we get an 'optimal' truncation order n op t = 3, 
with a 10% error on the exact value of the critical energy. We may wonder 
if there is a way to speed up the convergence rate: this can be done with a 
resummation method like the continued fraction [26] , It can be shown that, 
for all values of 5 up to 0.3, n = 6 is enough to reach a precision comparable to 
the numerical error. For 6 — 0.4 we get an optimal truncation order n opt = 5, 
with a 3% error on the exact value of the critical energy. 

3.2 Periodic orbits in general position and boxlets 

In addition to the normal modes, each resonant normal form of the type 
(|25p admits a double family of resonant periodic orbits in general position 
usually called boxlets in galactic dynamics [SB] . They can be easily identified 
using the fact that the two 'angles' have a fixed phase relation given either 
by ip = or ip = ±ir. In addition to the already mentioned 1:2 resonant 
banana, we have the 1:1 'loop', the 2:3 'fish' and so forth [28]. Each of them, 
if stable, is surrounded by a family of quasi-periodic orbits usually inheriting 
the same nickname. In a given potential in the class of ([T]), several boxlets 
can be present at the same time, whereas each resonant normal form is able 
to correctly render only one type: even with this limitation, a knowledge of 
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the corresponding family is very useful. In particular, we can analytically 
compute the phase-space fraction occupied by the given family, an important 
information in the process of constructing self-consistent models. We illustrate 
the idea in the case of the loop family: the principle is the same for higher 
resonating boxlets but the computations much more involved. 

Usually the loop bifurcates from the minor-axis periodic orbit at energy 
lower than that of the banana bifurcation fron the major axis 16 : there is a 
regime in which the loop family and the boxes around the stable major-axis 
orbit coexist. To identify loops, we impose the condition that the Hamiltonian 
flow generated by the 1:1 version of the reduced normal form K(lZ,ip;£,5) 
has a fixed point in TZ = TZl , ip = tt. Using relations ([26]) and |28|) and the 
value of the detuning S = q — 1, this solution fixes the actions on the closed 
loop: JiL(£,q) and J 2 l = £ — J\h{£,q)- On the periodic orbit, it is possible 
to find a relation between £ and the true energy £ in a form analogous to 
the expansion ([50]) . We can then express the actions as a series in E and, 
exploiting their geometric meaning, produce an estimate of the fraction of 
phase space occupied by the loops and the boxes. Truncating the series at 
first order, in the logarithmic case the results are [UJ 

_ J 1L (E, g) _ 2(-3 + 3g - 5g 2 + 5g 3 ) + E(9 -9g+ llg 2 - 3g 3 ) 
tLoop- s ^ - i3 _ 2 q + 3q 2 )(-2(l-q) 2 + E{3-2q + 3q 2 )) 

(34) 

and Jbox = 1 — fhoop- These predictions and the corresponding one for the 
banana family [17] agree very well with numerical estimates [28j up to energy 
values much greater than that corresponding to the harmonic core of the 
potential. 

3.3 Singular limits 

The potentials considered up to now are assumed to be analytic in the origin. 
However, we know that realistic models should include a singularity related to 
a density 'cusp' and/or a central point mass [30] • In the examples of the form 
([3l4p is implicitly assumed the use of adimensional coordinates by introducing 
a 'core radius' R c which can be put equal to 1 without less of generality. In 
the limit R c — *■ 0, those examples reduce to members of the familiy of singular 
scale-free potentials [31] 

V a (s)=As a , -l/2<a<l, Aa>0. (35) 

The singular conical potential is given by a = 1/2 and A = 1 while the 
singular logarithmic potential corresponds to the limit a — > with A = 1/2. 

It is tempting to try to extract information concerning the scale-free sin- 
gular limit from our analytical setting based on series expansions. Formally, 
this operation should be hindered by the lack of a series representation of the 
singular potential. However, we may nonetheless 'force' our approximate inte- 
grals of motion to play their role in the singular limit too and try our chance 
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by constructing a Poincare surface of section by using the approximate inte- 
gral I(x, y,p xi p y ; q) given by (f!ZtJ|) . In view of the scale invariance, we fix the 
energy level E and construct, e.g., a y-p y surface of section by means of the 
intersection of the function I{x 1 y,p y ] Eo, q) with the x — hyperplane. The 
level curves of the function 

F(y, Py )=I(0,y,p v ;E ,q) (36) 

give the invariant curves on the section. In the singular logarithmic case with 
E = and q = 0.7 we get, quite surprisingly, acceptable results [T7j. In the 
section constructed by using the approximate integral I^ la > related to the 1:1 
resonance and obtained by truncating the series (|20p at order 6, it gives the 
family of loops around the stable periodic orbit at y ~ 0.56 in good agreement 
with numerical data [28) . A similar result is given with the section obtained 
by using the approximate integral /( 1:2 ) again truncated at order 6: it gives 
the family around the stable banana at y ~ 0.16 and boxes around it. 

4 Comments and outlook 

As any analytical approach, this method has the virtue of embodying in (more 
or less) compact formulas simple rules to compute specific quantities, giving 
a general overview of the behavior of the system. In the case in which a non- 
integrable system has a regular behavior in large part of its phase space, a very 
conservative strategy, like that of truncating at a low order including the res- 
onance, provides sufficient qualitative and quantitative agreement with other 
more accurate but less general approaches. In our view, the most relevant 
limitation, common to all perturbation methods, is due to the intrinsic struc- 
ture of a single-resonance normal form. However, we remark that the regular 
dynamics of a non-integrable system can be imagined as a superposition of 
very weakly interacting resonances. If we are not interested in thin stochastic 
layers, each portion of phase space associated with a given resonance has a 
fairly good alias in the corresponding normal form. 

There are several lines of developement of this line of research; we men- 
tion a few of them: to extend the asymptotic analysis of perturbation series 
representing the building blocks of phase space (actions, frequencies, etc.); to 
devise suitable coordinate transformations to enable the investigation of cuspy 
potentials and/or central 'black holes'; to apply the normalization algorithm 
to three degrees of freedom systems with and without rotation. 
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